﻿// OCLMatrixMultiply.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "OpenCLKernel.h"
#include <math.h>
#include <time.h>


// multiply.cl
const char* _kernel_multiply_01_ = KERNEL(
    __kernel void multiply_01(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    const int i = get_global_id(0), j = get_global_id(1);
    int k;
    float sum = 0.0f;
    for (k = 0; k < colA; k++)
        sum += inputA[i * colA + k] * inputB[k * colB + j];
    outputC[i * colB + j] = sum;
    return;
}
);

const char* _kernel_multiply_02_ = KERNEL(
    __kernel void multiply_02(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    float Asub[1024];
    __local float Bsub[1024];
    const int i = get_global_id(0);
    const int li = get_local_id(0);
    const int nloc = get_local_size(0);
    int k, j;
    float sum;
    for (k = 0; k < colA; k++)
        Asub[k] = inputA[i * colA + k];
    for (j = 0; j < colB; j++)
    {
        for (k = li; k < colA; k = k + nloc)
            Bsub[k] = inputB[k * colB + j];
    barrier(CLK_LOCAL_MEM_FENCE);
        sum = 0.0f;
        for (k = 0; k < colA; k++)
            sum += Asub[k] * Bsub[k];
        outputC[i * colB + j] = sum;
    }
    return;
}
);

const char* _kernel_multiply_03_ = KERNEL(
    __kernel void multiply_03(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    __local float Asub[16][16], Bsub[16][16];
    const int li0 = get_local_id(0), li1 = get_local_id(1);
    const int gr0 = get_group_id(0), gr1 = get_group_id(1);
    const int nTiles = colA / 16;
    int gk, k;
    float acc = 0.0f;

    for (gk = 0; gk < nTiles; gk++)
    {
        Asub[li0][li1] = inputA[(gr0 * 16 + li0) * colA + gk * 16 + li1];
        Bsub[li1][li0] = inputB[(gk * 16 + li0) * colB + gr1 * 16 + li1];
        barrier(CLK_LOCAL_MEM_FENCE);

        for (k = 0; k < 16; k++)
            acc += Asub[li0][k] * Bsub[li1][k];
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    outputC[(gr0 * 16 + li0) * colB + gr1 * 16 + li1] = acc;
    return;
}
);

const char* _kernel_multiply_04_ = KERNEL(
    __kernel void multiply_04(__global float* inputA, __global float* inputB, __global float* outputC, int colA, int colB)
{
    __local float4 Asub[16][16], Bsub[16][16];
    const int li0 = get_local_id(0), li1 = get_local_id(1);
    const int gr0 = get_group_id(0), gr1 = get_group_id(1);
    const int nTiles = colA / 16;
    int gk, k, j;
    float4 acc[4] = { (float4)0.0f };
    
    for (gk = 0; gk < nTiles; gk++)
    {
        Asub[li0][li1] = (float4)(inputA[(gr0 * 16 * 4 + li0) * colA + gk * 16 + li1],
            inputA[(gr0 * 16 * 4 + 16 + li0) * colA + gk * 16 + li1],
            inputA[(gr0 * 16 * 4 + 16 * 2 + li0) * colA + gk * 16 + li1],
            inputA[(gr0 * 16 * 4 + 16 * 3 + li0) * colA + gk * 16 + li1]);
        Bsub[li1][li0] = vload4(((gk * 16 + li0) * colB + gr1 * 16 * 4 + li1 * 4) / 4, inputB);
        barrier(CLK_LOCAL_MEM_FENCE);

        for (k = 0; k < 16; k++)
        {
            acc[0] += Asub[li0][k].x * Bsub[li1][k];
            acc[1] += Asub[li0][k].y * Bsub[li1][k];
            acc[2] += Asub[li0][k].z * Bsub[li1][k];
            acc[3] += Asub[li0][k].w * Bsub[li1][k];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    
    for (j = 0; j < 4; j++)
        vstore4(acc[j], ((gr0 * 16 * 4 + 16 * j + li0) * colB + gr1 * 16 * 4 + li1 * 4) / 4, outputC);
    return;
}
);


// 4096 x 1024 * 1024 x 2048，暂不考虑除不尽等情况，可补0简单解决
const int rowA = 4096, colA = 1024, colB = 2048;


bool floatEq(const float a, const float b)// 相等返回 1
{
    if (b == 0)
        return fabs(a) < 0.001;
    return fabs(1.0 * a / b - 1) < 0.001;
}


int main()
{
    float* A = (float*)malloc(sizeof(float) * rowA * colA);
    float* B = (float*)malloc(sizeof(float) * colA * colB);
    float* C = (float*)malloc(sizeof(float) * rowA * colB);
    srand(2);
    for (int i = 0; i < rowA * colA; A[i] = rand() & 0xF, i++);
    for (int i = 0; i < colA * colB; B[i] = rand() & 0xF, i++);
    for (int i = 0; i < rowA * colB; C[i] = 0, i++);

    cl_int status;

    OpenCLKernel* m_pOCL = new OpenCLKernel;
    if (0 != m_pOCL->OpenCLInitialize())
    {
        delete m_pOCL;
        m_pOCL = nullptr;
        return EXIT_FAILURE;
    }

    cl_command_queue queue = clCreateCommandQueue(m_pOCL->GetContext(), m_pOCL->GetDevices()[0], NULL, &status);
    GET_CL_ERROR(status, m_pOCL);

    //编译
    size_t nSourceSize = NULL;
    cl_program program = NULL;
    cl_kernel kernel = NULL;
    printf("Choose the way of multiply.\n");
    int sel = 0;
    scanf_s("%d", &sel, 1);
    switch (sel)
    {
    case 1:
        nSourceSize = { strlen(_kernel_multiply_01_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_01_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        kernel = clCreateKernel(program, "multiply_01", &status);
        GET_CL_ERROR(status, m_pOCL);
        break;
    case 2:
        nSourceSize = { strlen(_kernel_multiply_02_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_02_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        kernel = clCreateKernel(program, "multiply_02", &status);
        GET_CL_ERROR(status, m_pOCL);
        break;
    case 3:
        nSourceSize = { strlen(_kernel_multiply_03_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_03_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        kernel = clCreateKernel(program, "multiply_03", &status);
        GET_CL_ERROR(status, m_pOCL);
        break;
    case 4:
        nSourceSize = { strlen(_kernel_multiply_04_) };
        program = clCreateProgramWithSource(m_pOCL->GetContext(), 1, &_kernel_multiply_04_, &nSourceSize, &status);
        GET_CL_ERROR(status, m_pOCL);
        status = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
        GET_CL_ERROR(status, m_pOCL);
        kernel = clCreateKernel(program, "multiply_04", &status);
        GET_CL_ERROR(status, m_pOCL);
        break;
    case 0:
    case -1:
    default:
        return EXIT_FAILURE;
    }

    //读内存
    cl_mem bufferA = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * rowA * colA, A, &status);
    GET_CL_ERROR(status, m_pOCL);
    cl_mem bufferB = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(float) * colA * colB, B, &status);
    GET_CL_ERROR(status, m_pOCL);
    cl_mem bufferC = clCreateBuffer(m_pOCL->GetContext(), CL_MEM_WRITE_ONLY, sizeof(float) * rowA * colB, NULL, &status);
    GET_CL_ERROR(status, m_pOCL);

    clock_t time_start = clock();

    //参数赋值
    status = clSetKernelArg(kernel, 0, sizeof(cl_mem), &bufferA);
    GET_CL_ERROR(status, m_pOCL);
    status = clSetKernelArg(kernel, 1, sizeof(cl_mem), &bufferB);
    GET_CL_ERROR(status, m_pOCL);
    status = clSetKernelArg(kernel, 2, sizeof(cl_mem), &bufferC);
    GET_CL_ERROR(status, m_pOCL);
    status = clSetKernelArg(kernel, 3, sizeof(int), &colA);
    GET_CL_ERROR(status, m_pOCL);
    status = clSetKernelArg(kernel, 4, sizeof(int), &colB);
    GET_CL_ERROR(status, m_pOCL);

    // 注意不同函数需要调整工作组网格参数
    size_t * globalSize = NULL;
    size_t * localSize = NULL;
    cl_uint work_dim = 0;
    switch (sel)
    {
    case 1:
    case 3:
        work_dim = 2;
        globalSize = new size_t[2];
        globalSize[0] = rowA;
        globalSize[1] = colB;
        localSize = new size_t[2];
        localSize[0] = TILE_DIM;
        localSize[1] = TILE_DIM;
        break;
    case 2:
        work_dim = 1;
        globalSize = new size_t[1];
        localSize = new size_t[1];
        globalSize[0] = rowA;
        localSize[0] = 256;
        break;
    case 4: //此算法的size皆以float4为单位
        work_dim = 2;
        globalSize = new size_t[2];
        globalSize[0] = rowA / 4;
        globalSize[1] = colB / 4;
        localSize = new size_t[2];
        localSize[0] = TILE_DIM;
        localSize[1] = TILE_DIM;
        break;
    default:
        return EXIT_FAILURE;
    }
    
    //执行
    status = clEnqueueNDRangeKernel(queue, kernel, work_dim, NULL, globalSize, localSize, 0, NULL, NULL);
    GET_CL_ERROR(status, m_pOCL);
    status = clFinish(queue);
    GET_CL_ERROR(status, m_pOCL);

    clock_t time_end = clock();
    printf("\nTime kernel : %d ms\n", time_end - time_start);

    //写内存
    status = clEnqueueReadBuffer(queue, bufferC, CL_TRUE, 0, sizeof(float) * rowA * colB, C, 0, NULL, NULL);
    GET_CL_ERROR(status, m_pOCL);

    // 返回并检查结果
    printf("Checking...\n");
    int correct = 1;
    for (int i = 0; i < rowA && correct; i++)
    {
        for (int j = 0; j < colB && correct; j++)
        {
            float tempSum = 0.0f;
            for (int k = 0; k < colA; k++)
            {
                tempSum += A[i * colA + k] * B[k * colB + j];
            }
            if (!floatEq(tempSum, C[i * colB + j]))
            {
                printf("Error at [%d, %d], calculation: %f, reference: %f\n", i, j, C[i * colB + j], tempSum);
                correct = 0;
            }
        }
    }
    if (correct)
        printf("Result correct.\n");

    //结束并释放
    if (globalSize != NULL)
        delete[] globalSize;
    if (localSize != NULL)
        delete[] localSize;
    clReleaseMemObject(bufferA);
    clReleaseMemObject(bufferB);
    clReleaseMemObject(bufferC);
    clReleaseKernel(kernel);
    clReleaseProgram(program);
    clReleaseCommandQueue(queue);
    delete m_pOCL;
    free(A);
    free(B);
    free(C);
    return 0;
}
